library(Boruta)
library(ggplot2)
library(ggpubr)
library(pheatmap)
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
#expr <- read.table("easy_input_expr.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
#tumsam <- colnames(expr)[substr(colnames(expr),11,12) == "01"]
#expr <- expr[,tumsam]
#indata <- t(scale(t(log2(expr + 1)))) # z-score表达谱

# 加载亚型结果
#annCol <- read.table("easy_input_cluster.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
table(annCol$Clust2)

# 加载Figure201ClusterCorrelation的ICI signature gene结果
outTab <- read.table("ouput_ICIsignatureGene.txt",sep = "\t",row.names = NULL,header = T,stringsAsFactors = F,check.names = F)
table(outTab$direct) #


set.seed(20201024)
dat.boruta <- as.data.frame(t(indata[outTab$gene,rownames(annCol)]))
borutafit <- Boruta(x = as.matrix(dat.boruta), 
                    y = as.factor(annCol$`Gene cluster`), # multiclassification
                    doTrace = 2,
                    maxRuns = 100,
                    ntree = 500)
boruta_fea <- attStats(borutafit)
boruta_fea <- rownames(boruta_fea[which(boruta_fea$decision == "Confirmed"),])
boruta.all <- outTab[which(outTab$gene %in% boruta_fea),]
table(boruta.all$direct)

# 拆分降维后的AB签名
boruta.A <- boruta.all[boruta.all$direct == "A",]
boruta.B <- boruta.all[boruta.all$direct == "B",]




###PCA_score
expr.A <- indata[boruta.A$gene,rownames(annCol)]
pca.A <- prcomp(t(expr.A), scale = F, center = F) # 如果数据没有标准化，这里都要设置为TRUE
pca1.A <- pca.A$x[,1] # 取出第一主成分

expr.B <- indata[boruta.B$gene,rownames(annCol)]
pca.B <- prcomp(t(expr.B), scale = F, center = F) #如果数据没有标准化，这里都要设置为TRUE
pca1.B <- pca.B$x[,1] # 取出第一主成分

ICI.score <-abs( pca1.A - pca1.B )# 主成份相减得到ICI得分
ICI.outtab <- data.frame(samID = rownames(annCol),
                         pca1.A = pca1.A[rownames(annCol)],
                         pca1.B = pca1.B[rownames(annCol)],
                         ICI.score = ICI.score[rownames(annCol)],
                         subtype = annCol$`Gene cluster`,
                         stringsAsFactors = F)
# 输出到文件
write.csv(ICI.outtab,"output_ICI_score.csv")


# 画图对比ABC三类的ICI score

参考FigureYa162boxViolin的画法


ICI.outtab<-merge(ICI.outtab,tcga_gsva[,c(1,10)])


ggplot(data = ICI.outtab,aes(x = subtype, y = ICI.score, fill = subtype))+
  scale_fill_manual(values = c("#008ECB", "#EA921D", "#D14039")) + 
  geom_violin(alpha=0.4, position = position_dodge(width = .75),
              size=0.8, color="black") + # 边框线黑色
  geom_boxplot(notch = TRUE, outlier.size = -1, 
               color="black", lwd=0.8, alpha = 0.7)+ # 背景色透明化
  geom_point(shape = 21, size=2, 
             position = position_jitterdodge(), 
             color="black", alpha=1)+ # 边框线黑色
  theme_classic() +
  ylab(expression("ICI score")) +
  xlab("Subtype")  +
  theme(axis.ticks = element_line(size=0.2,color="black"),
        axis.ticks.length = unit(0.2,"cm"),
        legend.position = "none",
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10)) +
  stat_compare_means(method = "kruskal.test", label.y = max(ICI.score))
ggsave(filename = "ICI_score.pdf", width = 5, height = 5)
```

# Session Info
##ALL
clin<-d
colnames(ICI.outtab) [1]<-"ID"
clin$ID<-rownames(clin)
all_clin<-merge(clin,ICI.outtab[,c(1,4)],by="ID")




library(survminer)
library(survival)
library(survminer)
res.cut <- surv_cutpoint(all_clin, time = "futime",
                         event = "fustat",
                         variables = names(all_clin)[32],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
all_clin$group<-res.cat$ICI.score
write.csv(all_clin,"all_clin.csv")

Sur <- Surv(as.numeric(all_clin[,10]) ,as.numeric(all_clin[,11]))
sfit <- survfit(Sur ~ all_clin[,33],data=as.data.frame( all_clin[,33]) )

ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_score','Low_score'))+
  labs(x = "Days")

fitd <- survdiff(Surv(futime, fustat) ~ group,
                 data = all_clin,
                 na.action = na.exclude)
p.val <- 1-pchisq(fitd$chisq, length(fitd$n) - 1)


HR = (fitd$obs[1]/fitd$exp[1])/(fitd$obs[2]/fitd$exp[2])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))

metadata1<-data.frame(row.names = "TCGA",
                      HR=HR,
                      up95=up95,
                      low95=low95)




###GSE27831
GSE27831<-read.csv("GSE27831_select2.csv",header = T,row.names = 1,check.names = F)
GSE27831$ICI.score<-abs(GSE27831$UFC1-GSE27831$JUP)

res.cut <- surv_cutpoint(GSE27831, time = "days",
                         event = "event",
                         variables = names(GSE27831)[5],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
GSE27831$group<-res.cat$ICI.score


fitd <- survdiff(Surv(days, event) ~ group,
                 data = GSE27831,
                 na.action = na.exclude)
p.val <- 1-pchisq(fitd$chisq, length(fitd$n) - 1)


HR = (fitd$obs[1]/fitd$exp[1])/(fitd$obs[2]/fitd$exp[2])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))

metadata2<-data.frame(row.names = "GSE27831",
                      HR=HR,
                      up95=up95,
                      low95=low95)



Sur <- Surv(as.numeric(GSE27831[,1]) ,as.numeric(GSE27831[,2]))
sfit <- survfit(Sur ~ GSE27831[,6],data=as.data.frame( GSE27831[,6]) )

ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_score','Low_score'))+
  labs(x = "Days")


###GSE22138
GSE22138<-read.csv("GSE22138_select2.csv",header = T,row.names = 1,check.names = F)
GSE22138$ICI.score<-abs(GSE22138$UFC1-GSE22138$JUP)
res.cut <- surv_cutpoint(GSE22138, time = "days",
                         event = "event",
                         variables = names(GSE22138)[5],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
GSE22138$group<-res.cat$ICI.score


fitd <- survdiff(Surv(days, event) ~ group,
                 data = GSE22138,
                 na.action = na.exclude)
p.val <- 1-pchisq(fitd$chisq, length(fitd$n) - 1)


HR = (fitd$obs[1]/fitd$exp[1])/(fitd$obs[2]/fitd$exp[2])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))

metadata3<-data.frame(row.names = "GSE22138",
                      HR=HR,
                      up95=up95,
                      low95=low95)



Sur <- Surv(as.numeric(GSE22138[,1]) ,as.numeric(GSE22138[,2]))
sfit <- survfit(Sur ~ GSE22138[,6],data=as.data.frame( GSE22138[,6]) )

ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_score','Low_score'))+
  labs(x = "Days")

##

###GSE84976
GSE84976<-read.csv("GSE84976_select2.csv",header = T,row.names = 1,check.names = F)
GSE84976$ICI.score<-abs(GSE84976$UFC1-GSE84976$JUP)
res.cut <- surv_cutpoint(GSE84976, time = "days",
                         event = "event",
                         variables = names(GSE84976)[5],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
GSE84976$group<-res.cat$ICI.score


fitd <- survdiff(Surv(days, event) ~ group,
                 data = GSE84976,
                 na.action = na.exclude)
p.val <- 1-pchisq(fitd$chisq, length(fitd$n) - 1)


HR = (fitd$obs[1]/fitd$exp[1])/(fitd$obs[2]/fitd$exp[2])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))

metadata4<-data.frame(row.names = "GSE84976",
                      HR=HR,
                      up95=up95,
                      low95=low95)



Sur <- Surv(as.numeric(GSE84976[,1]) ,as.numeric(GSE84976[,2]))
sfit <- survfit(Sur ~ GSE84976[,6],data=as.data.frame( GSE84976[,6]) )

ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_score','Low_score'))+
  labs(x = "Days")


###GSE44295
GSE44295<-read.csv("GSE44295_select2.csv",header = T,row.names = 1,check.names = F)
GSE44295$ICI.score<-abs(GSE44295$UFC1-GSE44295$JUP)
res.cut <- surv_cutpoint(GSE44295, time = "days",
                         event = "event",
                         variables = names(GSE44295)[5],
                         minprop = 0.3)

res.cat <- surv_categorize(res.cut)
GSE44295$group<-res.cat$ICI.score


fitd <- survdiff(Surv(days, event) ~ group,
                 data = GSE44295,
                 na.action = na.exclude)
p.val <- 1-pchisq(fitd$chisq, length(fitd$n) - 1)


HR = (fitd$obs[1]/fitd$exp[1])/(fitd$obs[2]/fitd$exp[2])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))

metadata5<-data.frame(row.names = "GSE44295",
                      HR=HR,
                      up95=up95,
                      low95=low95)



Sur <- Surv(as.numeric(GSE44295[,1]) ,as.numeric(GSE44295[,2]))
sfit <- survfit(Sur ~ GSE44295[,6],data=as.data.frame( GSE44295[,6]) )

ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_score','Low_score'))+
  labs(x = "Days")


###E-MTAB-4097
E_MTAB<-read.csv("E-MTAB-4097.csv",header = T,row.names = 1,check.names = F)
E_MTAB$ICI.score<-abs(E_MTAB$UFC1-E_MTAB$JUP)

res.cut <- surv_cutpoint(E_MTAB, time = "days",
                         event = "event",
                         variables = names(E_MTAB)[5],
                         minprop = 0.2)

res.cat <- surv_categorize(res.cut)
E_MTAB$group<-res.cat$ICI.score


fitd <- survdiff(Surv(days, event) ~ group,
                 data = E_MTAB,
                 na.action = na.exclude)
p.val <- 1-pchisq(fitd$chisq, length(fitd$n) - 1)


HR = (fitd$obs[1]/fitd$exp[1])/(fitd$obs[2]/fitd$exp[2])
up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))
low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/fitd$exp[1]+1/fitd$exp[2]))

metadata6<-data.frame(row.names = "E-MTAB-4097",
                      HR=HR,
                      up95=up95,
                      low95=low95)



Sur <- Surv(as.numeric(E_MTAB[,1]) ,as.numeric(E_MTAB[,2]))
sfit <- survfit(Sur ~ E_MTAB[,6],data=as.data.frame( E_MTAB[,6]) )

ggsurvplot(sfit, 
           conf.int=F, #置信区间
           #fun="pct",
           pval=TRUE,
           palette = c("#EA921D", "#008ECB"),
           pval.method = T,
           risk.table =T, 
           ncensor.plot = F,
           surv.median.line="hv",
           legend.labs=c('High_score','Low_score'))+
  labs(x = "Days")
